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NUMERICAL CALCULATION OF BOUND ARY -INDUCED 
INTERFERENCE IN SLOTTED OR PERFORATED WIND TUNNELS 
INCLUDING VISCOUS EFFECTS IN SLOTS 

By James D. Keller 
Langley Research Center 

SUMMARY 

A numerical method is presented for calculating the incompressible boundary - 
induced interference in wind tunnels of rectangular cross section with slotted or perfo- 
rated walls. The method includes a wall representation which is capable of satisfying a 
generalized homogeneous boundary condition including the effects of viscosity within the 
slots. The effects of viscosity in the slots are found to be very significant. The method 
allows for a variation in the boundary conditions along the tunnel walls. The model can 
be any configuration and can be located anywhere in the test section. The interference 
can be computed at any point in the test section. 

INTRODUCTION 

In order to obtain accurate wind-tunnel data, the measured quantities must often be 
corrected to account for the interference caused by the wind-tunnel boundaries. Theo- 
retical methods are presently available for predicting the interference due to the tunnel 
walls in certain cases. The analytical methods are limited to infinite -length test sections 
with constant wall properties in the tunnel stream direction. Some methods are limited 
as to model size, position, and load distribution. A numerical method for calculating 
the boundary -induced interference in ventilated wind tunnels is presented in reference 1. 
The method consists of dividing the tunnel walls into rectangular elements which are each 
represented by a source distribution. A matrix equation is then solved to find the source 
strengths which allow the boundary conditions to be satisfied at the centroid of each ele- . 
ment. In reference 1 each element was represented by a source distribution of constant 
strength over the element. This representation is particularly well suited to satisfying 
an ideal slotted-wall boundary condition. The ideal slotted-wall boundary condition, how- 
ever, is only a special case of a more general boundary condition which can include the 
effects of viscosity within the slots. 

The present investigation deals specifically with a modified representation for the 
tunnel walls which is suitable for the satisfaction of the more general boundary condition 


including the effects of viscosity in the slots. The method presented is limited to incom- 
pressible flow and cannot handle the usual assumption of a test section which extends to 
infinity upstream and downstream of the model. The method also requires the experi- 
mental determination of one of the parameters in the boundary condition. The method has 
broad applicability, however, because the boundary conditions on the tunnel walls may 
vary almost without limit. The model representation is also quite general. The model 
may be located anywhere in the test section and at any orientation. A sample computer 
program used in making the calculations is given in an appendix. 

SYMBOLS 

a effect of one element on another 

b effect of model on an element 

c l’ c 2 ,c 3 ,c 4 coefficients in equation (4) 
d distance between slot centers 

1 slot parameter 

N total number of elements 

n direction normal to wall 

R restriction parameter 

s wing span 

t slot width 

w w upwash velocity caused by tunnel walls 

x,y,z Cartesian coordinates 

r m circulation of model 

6 lift interference factor 

|,77,£ Cartesian coordinates 

2 



a source distribution strength 



q> perturbation velocity potential function 

cp* velocity potential function for an element divided by o’ for the element 

Subscripts: 

i at ith element 

j at jth element 

L downstream end of each source distribution 

m due to model 

w due to tunnel walls 


ANALYSIS 

General Statement of Problem 

The governing equation used in the analysis of incompressible wind-tunnel inter- 
ference is 

i?£ + J?i£ + i!£ = o (1) 

ax 2 ay 2 dz 2 

where cp is the perturbation velocity potential function for the entire flow field. Let 
cp = <p m + (p^ where <p m is the potential function of the disturbances due to the model 
in free air and <p w is the potential function of the additional flow due to the tunnel walls. 
If cp m is taken as a known solution of equation (1) which approximates the flow field at 
a distance from the model in free air, then <p w can be determined by the fact that cp 
must satisfy certain boundary conditions at the tunnel walls. The objective in determining 
cp w is to be able to calculate the change in the free -stream conditions caused by the tun- 
nel walls. Since cp m needs to be known only on the tunnel walls, any inaccuracies in 
cp m near the model will have no effect on the determination of <p w . 
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Boundary Conditions 


The boundary condition to be satisfied at a solid boundary is that there can be no 
flow through the boundary; that is 

0 

9n 

where n is the direction normal to the wall (positive outwards). The boundary condi- 
tion to be satisfied at an open jet boundary is that there be no pressure difference across 
the boundary. This boundary condition can be approximated by (ref. 2) 

o 

ax 

Reference 3 gives the homogeneous boundary condition to be satisfied at a perforated wall 
as 

l£ + ii£= o 

ax R an 

where R is a restriction parameter which relates the pressure difference across the 
wall to the flow through the wall. In practice, R must be determined experimentally 
for a given wall. A detailed discussion of the restriction parameter for both porous and 
perforated walls can be found in reference 4. 

The homogeneous boundary condition to be satisfied at a wall with several longitu- 
dinal slots is given in reference 5 as 

cp + l = 0 (2) 

where l is a slot parameter given by 


l dincsc( 2d ) 

where t is the slot width and d is the distance between slot centers. This slot param- 
eter was derived on the basis of two-dimensional flow, and it is assumed that it can be 
applied at each location in the tunnel even if the slot width varies. Equation (2) can be 
differentiated with respect to x to give 



For constant slot width, equation (3) becomes 


| £ + ; 

ax 


d ^Cp _ 


dx 9n 


which is the form given by many authors. 

The ideal slotted-wall boundary condition was derived on the basis of inviscid flow 
at the slots. In reference 6 the addition of another term to account for the effects of vis- 
cosity in the slots is suggested. The boundary condition is then 
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ie + lig + z _ 9 ‘ J£- = o 

dx R 9n 3x 9n 


If the coefficient l/R is replaced by i- + this boundary condition will also apply to 

Xv C/X 

walls which do not have constant slot width. As R approaches infinity, this boundary 
condition approaches that for an ideal slotted wall. However, as pointed out in refer- 
ence 7, the ideal condition is not always valid. Experiments have shown that for typical 
test -section configurations, R can be of the order of unity (ref. 8). This value of R 
will have a very pronounced effect on the interference in the test section. Thus, it is 
important to consider the effects of viscosity in the slots and retain the additional term 
in the boundary condition. Note that no theoretical method exists for determining the 
value of R for a particular test section. It must be determined experimentally. This 
experimental determination of the restriction parameter might be quite difficult because 
for a tunnel with varying slot width the restriction parameter may vary with position also. 


In this paper, a general boundary condition of the form 


c X <P 


+ c.^ + c.ise + c. T 2-¥-=o 


9x 


3 9n 4 ^ 9 n 


( 4 ) 


is considered. This boundary condition contains all previous conditions as special cases 
as shown in the following table: 


Type of boundary 
condition 

C 1 

c 2 

c 3 

c 4 

Closed wall 

0 

0 

1 

0 

Open jet 

0 

1 

0 

0 

Perforated wall 

0 

1 

1 

R 

0 

Ideal slotted wall 

1 

0 

l 

0 

(integrated form) 



se 

dx. 


Ideal slotted wall 

0 

1 

l 

(differentiated form) 




Slotted wall including 

0 

1 

— + — 

i ax r 

i 

l 

viscosity in slots 





Additional discussions of these boundary conditions may be found in references 4, 6, and 7. 

Representation of Tunnel Walls 

In order to satisfy the homogeneous boundary condition, the tunnel walls are divided 
into longitudinal strips and each strip is divided into a number of rectangular elements. 
The boundary condition will be satisfied at the centroid of each element. The coordinate 
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system to be used has the X-axis extending along the tunnel center line, with the positive 
direction being the tunnel stream direction. The Z-axis is positive upwards, and the 
Y-axis is chosen so that the coordinate system is a right-handed system. In reference 1 
each tunnel -wall element was represented by a source distribution of constant strength 
over the element. This representation is particularly suited to the satisfaction of the 
ideal slotted-wall boundary condition in integrated form (eq. (2)) because for this case 
the matrix of influence coefficients is diagonally dominant. This diagonal dominance also 
holds for the perforated-wall boundary condition for small values of R. However, for 
the ideal slotted-wall boundary condition in differentiated form or for large values of R 
in the general slotted-wall and perforated -wall boundary conditions, this representation 
would lead to elements on the diagonal of the matrix of influence coefficients which are 
either zero or very small. This nearly singular matrix leads to numerical difficulties 
and inaccuracies in trying to solve the resulting system of equations. In order to avoid 
these difficulties, let each tunnel -wall element be represented by a source distribution 
over the element and downstream of the element at least to the end of the strip. The 

j 

strength of the source distribution a varies linearly, with a slope o' = 22, on the ele- 
ment itself and then remains constant downstream of the element. This representation 
is used so as to have the source strength continuous along a strip and still have only one 
unknown (the source strength slope) for each element. Thus, the source strength is zero 
at the upstream end of the strip and varies in linear segments along the length of the 
strip. If cp* is the potential function for a particular element divided by the source 
strength slope o' for that element, then 
N 

< 5 > 

j=l 

where N is the total number of elements. 

Consider an element in the top or bottom wall with corners as shown in figure 1. 



(£| i 1 ?! ^2’^l 


Figure 1.- Schematic of an element in top or bottom wall. 
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The potential function at a point (x,y,z) due to this source distribution is 


cp* = _f* 2 p 2 (4 ~ ^ K l rv 2 fe - ^l) d7 ? 

*1 ^ \/(x - 4) 2 + (y - 77)2 + (z - ?1 ) 2 ^2 *1 \l(x - £) 2 + (y - 77 ) 2 + "(z - ^) 2 

This potential function and its derivatives must be evaluated to satisfy the boundary con- 
ditions. For convenience in writing the equations, let Xj = x - x 2 = x - £ 2 > 

Xl = x - y t = y - Vi, y 2 = y - t? 2> an< l z \ = z ~ ?]_• The required equations are then 



( 6 ) 
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In order to find the effect of an element in a side wall, y and z must be inter- 
changed in equations (6) to (11). By examining equations (10) and (11) for Zj — 0, the 
effect of an element at its own centroid is found to be 

an 2 

and 


ax an 

At points of the same strip but downstream of the element 

- 2 ’(*2 - * l ) 

and 



ax an 


At all other elements of the wall in which an element is located, its contributions to 
d<p* / 8n and B^cp*/dx an are zero. 


Computation of Source Strength Slopes 


In order to compute the source strength slopes required to satisfy the boundary con- 
ditions at the centroid of each element, a matrix equation is needed to express these 
boundary conditions. Let a^ be the effect at the centroid of the ith element due to the 

source distribution corresponding to the jth element fa = c^cp* + C2 + Cg 

g2 (7?*\ 

— 1 Let be the effect of the model at the centroid of the ith element 


+ c 4 ax an 

(b = c 1 (p m + c 


d( P m 


x '2 — + c 3 

the boundary condition is 

AS’ = -B 

where 


dcp. 


m 


a 2 cp 


an 


+ c 


4 dx Bn J ' ^'^ ien equation which expresses 


( 12 ) 


A [ a ij] 

M 

and 

b-H 

Equation (12) can be solved for the values of a! which can then be used to compute the 
interference potential due to the tunnel walls through the use of equation (5). 
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RESULTS AND DISCUSSION 


As a relatively simple example, consider the lift interference due to a small lifting 
wing mounted in the center of a square test section with solid side walls and four equally 
spaced slots in the top and bottom walls. Each tunnel wall is divided into four strips of 
equal width and each strip is divided into 10 elements by cutting planes at x = -1.00, 
-0.70, -0.45, -0.25, -0.10, 0.00, 0.10, 0.25, 0.45, 0.70, and 1.00. (For convenience, the 
tunnel width is taken to be unity.) 

The wing is represented by a horseshoe vortex of circulation T m . The span s 
of the horseshoe vortex is assumed to be so small that it becomes a vortex doublet start- 
ing at (0,0,0). The perturbation velocity potential function cp m at a point (x,y,z) due to 
this representation of the model is given by 


so that 


and 


^m 


r m s z 

47T y2 _j_ z 2 


1 + 


i/ 2 2 2 , 

Vx + y + z “j 


1 

a< ?m_ 

1 z 




Tm s 

ax 

4tt / o o 
+ y^ + 

W 2 



1 

d( Pm _ 

-1 

... 2x 3 y_z 

+ 3xy 3 z 

+ 3xyz 3 


r m s 

9y 

47r(y 2 + z 2 ) 2 

(x2 

+ y2 + z 

2 j3/ 2 


1 

a 2 (p 

3 

yz 



r s 
m 

ax ay 

4w / 2 
[x* + y 




1 

a< ?m _ 

1 

F 3 

v 2.z 2 + ^ 

Z 2 + xy 4 

3 2 

- x z - 

2 2 
xy z 

r s 
m 

az 

47r(y 2 + z 2 ) 2 


(* 2 

| 9 

+ y + 2 

,2)3/2 

1 

02 ®m 

1 x 2 + y^ 

- 2z 2 



TrvnS 

m 

ax az 

(x 2 + y 2 - 

2) 5 /2 
f z “•) 




These quantities are used on the right-hand side of equation (12) which is then solved for 
a jy/r m s - The values of a!y/r m s are suitable for the computation of the upwash velocity 
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1 


at any point in the test section by summing the velocity due to each 


w 


w 


r m s 


m 


element. 


9< ?w 
S dz 

The lift interference factor (ref. 5) is then 


6 = 4 


w 


w 


2 r s 
m 


Figure 2 shows the lift interference factor at the center of the tunnel as a function 
of the ratio of the slot width to the distance between slot centers. The results were com- 
puted by using the ideal slotted-wall boundary condition in both integrated and differenti- 
ated forms. The results computed by using the integrated form of the boundary condition 



io 1 l J I I I I i I I i 

' 0 .02 .04 .06 .08 .10 .12 .14 .16 .18 .20 

Ratio of slot width to distance between slot centers ,t/d 

Figure 2.- Lift interference factor at vanishingly small-span wing 
in square tunnel with four slots in top and bottom walls. 


(solid-line curve) are the same as those computed by using the wall representation of 
reference 1. The difference in the lift interference factor when the boundary condition is 
used in the two different forms is due to the fact that the differentiated form of the bound- 
ary condition does not satisfy the additional requirement that there be no disturbance due 
to the tunnel walls at an infinite distance upstream of the model. However, the addition 
of elements to the upstream end of the test section rapidly eliminates the difference. 

When the test section starts three tunnel widths upstream of the model, the difference is 
very nearly zero. The reason for this can be seen more clearly in figure 3 which shows 
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Longitudinol distance along strip 


Figure 3.- Variation of source strength along wall. 

a typical variation of the source distribution strength along one of the strips of the top 
wall. It can be seen that the distribution which corresponds to the differentiated form of 
the boundary condition would have the same value at the centroid of each element as the 
other distributions if the source strength was increased by a constant amount. This con- 
stant shift corresponds to the source strength which would have been built up by the 
point x = -1 on an infinitely long test section. The error caused by using the differ- 
entiated form of the boundary condition comes about not so much because the source dis- 
tribution on the far upstream portion of the test section is not included, but rather because 
the source strength which would have been built up on the far upstream portion of an 
infinitely long test section is not included all along the rest of the test section. Putting 
additional elements farther upstream of the model almost completely eliminates this 
error, but requires a larger matrix which takes up more computer storage and time. On 
the downstream end, the source distribution which extends downstream from each element 
can be extended beyond the end of the last element if desired. 

Figure 4 shows the lift interference factor at the center of the tunnel with the gen- 
eral slotted-wall boundary condition used for several values of the restriction param- 
eter R. This figure shows that the effects of viscosity in the slots can be very 
significant. 
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Figure 4.- Lift interference factor at vanishingly small- span wing 
in square tunnel with four slots in top and bottom walls for 
several values of R. 


The results presented here are for the simple case of a small -span wing mounted in 
the center of a tunnel with constant slot width and constant restriction parameter. The 
method, however, is applicable to more general problems. The model representation can 
be quite general, including the case of a large -span swept wing with nonuniform loading 
located anywhere in the test section. The slot width and restriction parameter may vary 
with position on the boundary. The results presented here are also for the case of a test 
section which extends one tunnel width upstream and downstream of the model. This 
arrangement of tunnel -wall elements was used for several reasons. First, the computer 
time required to invert a matrix of this size (N = 160) is not too large (about 1 minute on 
a Control Data 6600 computer system). Second, the results can be compared directly 
with those of reference 1 which used the same arrangement of elements but a different 
source distribution to represent each element. Third, use of this tunnel length clearly 
shows the error caused by using the differentiated form of the ideal slotted-wall boundary 
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condition and not including the portion of the boundary far upstream of the model. 

Greater accuracy could be obtained by adding elements to the test section, particularly 
farther upstream of the model, and by letting the constant strength source distribution 
which trails downstream of each element to extend beyond the end of the last element to 
a large distance downstream of the model. In equations (7) to (11), x^ can be allowed 
to approach minus infinity and simplify the equations somewhat. This is not possible 
in equation (6), so when the ideal slotted-wall boundary condition is used in integrated 
form, the test section must be terminated at a finite distance downstream of the model. 
Although the present development is oriented toward tunnels of rectangular cross section, 
it could be extended to other cross-sectional shapes. 

The sample computer program given in the appendix is not the one used to compute 
the results shown here. It has additional elements farther upstream and downstream of 
the model for greater accuracy. In the sample program the elements extend from 
x = -2.6 to x = 1.8, and the constant strength portion downstream of each element 
extends to = 10. This program requires about 240 000g storage locations and about 
5 minutes on a Control Data 6600 computer system. 

CONCLUDING REMARKS 

A numerical method for calculating the boundary -induced interference in slotted or 
perforated wind tunnels has been presented. The method includes a wall representation 
which is capable of satisfying a generalized boundary condition including the effects of 
viscosity within the slots of a slotted-wall tunnel. The method is limited to incompress- 
ible flow and requires the experimental determination of one of the coefficients in the 
boundary condition. It is also limited to finite -length test sections. 

The method presented has broad applicability and allows for the boundary condition 
to vary over the test -section walls. This feature should aid in the design of test sections 
which have nearly constant interference over the space occupied by the model. The 
model representation is also quite general, including the case of a large -span swept wing 
with nonuniform loading located anywhere in the test section. The interference can also 
be computed anywhere in the test section. The effects of viscosity in the slots are found 
to be very significant. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., June 30, 1972. 
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APPENDIX 


SAMPLE FORTRAN PROGRAM 

TMI s A P° EN D ' X CONTAINS A SAMPLE FP?TRAN P co GR 4*4 FOP COMPUTING TEE I 'FT 
• ^TE^FEPEMCE FACT CP IN A WIND Tl' NFL OF r FO T AMGLIL A A CROSS SECTION WITH 
SLOTTED OR PERFORATED WALLS. T' P P r G P A M WAS WCTTTEN FOR JSF OF' CDF 6000 

SE-' ,T ES COMPUTERS. IT IS INTENDED ONLY AS A SAMPl.*'. MODIFI C4TT ONS MUST PF 
MADE TO THE PPnQtiM T N ORDER TO COMPUTE DIFFERENT FASES. IF IT IS DtS 1 P ED 
TO REDUCE THE M AT® T X SIZE. THE PFODPAM GIVEN TN D f F £F FNC F L r DU I D PE 
‘MODIFIED USING T*E ARITHMETIC STATEMENT FlJNCTIf MS GIVEN ^ . FORMAL INPI'TS 
FAN RE FOUND DM I INES 203 AMD 269. 

PROGRAM A 37 71 ( I NOUT » OUTPUT, T APF5 = INPUT , TAPF 6 =PUT CUT ) 

D I MENS TUN X I ( 2 56 ) , E TA ( 256 ) . Z ET A ( 256 ) , X U ( 256 ) , X t 2 < 256 ) , 3 T Ai ( 25 6 ) ? 

* ,FTA2( 256) ,7ET A1 (256) , ZFTA2 ( 256 ), ' (256,256) , 8 ( 256) . S IONA (256) 0 

* , C 1(256), C. 2(256), C3(2 56), 04(2 56) A 

DIMENSION XA(1li),YA(L3),WT( ID) S 

SGN(X)=S T GN( 1.0, X) 6 

t (Xl,X2,Yl,Y2,Z)=Y2/2.*SOFT(X2"X2+Y2*Y2+Z*Z)-Yl/2.*S9RT(X2'X2+YL* 7 

* Y l + Z*Z )- Y2/2.* SOP T ( XI *X 1 + Y2* Y2+ Z* Z ) +Y i /?. S Of T ( X i * X . *Y . * Y . +Z* 7 ) 3 

* + ( Z*Z-X?*X2 )/? .*6LCD ( AOS ( ( Y2+SOPT ( X2*X? + Y2*Y2 + Z*Z ) ) / ( Yl + SCR T ( X2 J ) 

* X2+Yl*Yl+Z*Z> ) ) )-( Z*Z-X1*X1 )/?.*ALOG( A HS ( ( Y 2 + SD > T ( X I *X 1 *-Y2 * Y2 + Z * Z . .» 

+ > >/ ( Yl+SOP T ( XI *X1+Y1*YI +Z*Z ) ) ) ) —XI *Y2^ M nr, ( ' RS ( ( X2 + SD’ T ( X2*X2+Y2* 11 

Y2 + Z*Z ) ) / ( X. + S0PT( Xl*Xi + Y2*Y2 + Z*Z) ) ) )*X. a v.*al'IG( iiiS( ( X?FSQ f T (X 2 * i? 

* X2+Y1* Yl +Z*Z ) >/(Xl+SQRT(Xl*Xl+Yl*Yl + Z*Z) ) I )-XO*Xl*U np(A5S( (Y2 + 13 

* SDR T ( XL J 'XL+Y2*Y2 + Z*Z) )/<Yi. + S3R MXI *XL + Y' .^Y . +7*Z) ) ) )-X0*Y2< •* 1 0G( 

* RS ( (Xl +SCPT(Xl*XL+Y2*Y2«-Z*Z) ) / ( X2 ♦SQR T ( X2*X2 Y 2 *Y 2 + Z -Z ) )) )+XD*Yl 15 

* 15 ( ABS( ( X( + SOFT ( XL*X( +Yi*Y >Z*Z ) ) / ( X 2 + SDT. t ( x?*X? + Y ,«Y l + Z*Z ) ) ) ) i .6 

* -XL * A3 S ( Z)*(-ATAN( X?*Y2/ABS (Z )/SQR T ( X?*X2*Y2*Y2+Z*Z ) ) «-AT AM X1< ! Y2/ 17 

* APS(Z)/SDkT(Xl*X' l . + Y2*Y2+Z*7) )+ATAN(X2' s Yi/A I 3S(Z) / SQ y T(X2*X2+Yi*Yl+ 13 

* Z-*Z>)-ATAN(X1*YI/ARS(Z)/$QPT(X1*X1*Y1*Y1+Z*Z) ) )-XD*ARS( Z )*( -A TAM 19 

* X' *Y2 / A3S( Z ) /S0RT( XL*XL+Y?TY2<-Z*Z ) ) + AT an ( X2 *- Y? / AR S ( L ) / SORT ( X2 «X 2 + 2 1 

* Y2*Y2 + Z* Z ) ) +AT AM( XL*Y 1/ “ P-S ( Z ) /SO p T ( XL*XL+Y 1*Y1+Z*Z ) )-ATAM( X2* Yl/ 21 

* APS(Z)/SOrtT(X2*X2+Y'.*Yl + Z*Z) ) ) 22 

S( XI, X2, Yl, Y2,Z ) =Y2/2.*SQRT ( X 2* X2 +Y 2* Y? + Z* Z )-Yl/2.*SQR T ( X2--X2 + Y1* 23 

* Y1 + Z*Z)-Y 2 / 2 .*SDPT( X1*X’. + Y2* Y?+Z*7 )+Y L / 2 . * S 01 T ( yi *x . +Y . * Y- + /.* Z ) 24 

* + ( Z -Z-X2*X 2 1/2 . *A|. CD ( A SS ( ( Y2+$0 P T( X2*X2-t-Y.? 5 ’Y?+Z*7. ) ) / ( Y 1 + SOP T ( X2*" 25 

* X 2 +Yl*Y 1+7 *Z ) ) ) ) -( Z*Z-Xl*x 1 ) /?. *Al OG ( A 3S ( ( Y 2 + SCFT < X . *X *Y?»Y 2 + Z*Z 26 

* ) ) / ( Y l +SQR T ( X l*Xl+Y l*Yl+Z* Z ) ) ) ) — X 1 *Y 2* A LOT- ( A P S ( ( X2+SQR T( X2*X2+Y2* 27 

* Y?+Z*Z ) ) / ( Xl + SOR T( X1*X1+Y2 *Y2 + Z*7. ) ) ) ) + X 1 *Y 1 *A| OD ( A 3$ ( ( X2+S0R T ( X2 * 23 

* X?+YI*Y1 +1*1 ) ) / ( XI + SQRT ( X1*X1+Y1*Y 1 + Z*Z ) ) ) J-XD^XL^M DD( AOS ( ( Y2+ 2R 

* SORT( XI *XL +Y2*Y2 + Z*Z ) ) /( Yl +SO r T( XL*XL + Y'.*Y 1 +Z* Z ) ) I ) -XD*Y2* Al.CG( 3 ) 

* ADS ( ( XH-SQPT ( XL*XL+Y2*Y?+Z*Z ) ) / ( X2+SQF T ( X2*X2 fYZS'Y 2 +Z*Z ) ) ) ) + XD*Y 1 31 

* *Al 3G( ABS( ( XL + SOFT (XL*Xt+Yl*Y'. + Z*7))/< X2 + S 0 F T(X2*X?+Y.*YL+Z*Z)))) 32 

D°DX(X1,X2»Y1,Y2.Z) =Y 2 *AL n G ( A B S t ( Xl+SOPT( XI * XI + Y2*Y2 ^Z *Z ) )/ ( X2+ 3 3 

* SOR T ( X2*X2+Y2*Y2+Z*Z ) ) ) )+Y 1* 'LOG ( APS ( ( X2. + S0F T ( X?*X2+Y . *Y i.+ Z A Z ) )/ 34 

* ( XI +SORT (X 1*-X1+Y1*YH-Z*Z ) ) ) ) + X1*AL^F- ( ' BS ( ( Y1 + SQRT( K2*X?.*- Yl* Yl + 7*Z 35 

* ) )/(Y2 + SQRT(X2*X?+Y2*Y2+Z*Z) )*( Y2 + S0R T ( X l'-X 1+Y 2*Y 2-r Z* Z ) ) / ( Yl + SQF T 36 

* ( X1*X1 +Y 1*Y1+Z*Z ) ) ) ) +XD*AL no ( AB S ( ( Y 1 + SQR T ( XL ■i'XL+Yl’.” Y1 + Z*Z ) ) / ( Y2+ 37 

* SQRT( XL^XI +Y2*Y2 + Z*Z ) )*( Y2 + S0RT (X2*X2 + Y2*Y?+Z-Z ))/ ( Y L + SORT ( X2*X?«- 33 

Yl^Yl+Z+Z) ) ) ) 39 

* +ABS(Z)*(ATAM( X2*Y2/ABS(Z) / SDRT(X 2 *X 2 «-Y 2>! t Y2 +Z':-Z) ) - AT an ( x 1< Y2 / F 8 S < 4 ) 

* Z ) / SDR T ( X1*X1<-Y2*Y?+Z*Z ) )-ATAN( X2*YL/ABS( Z) /SQ»T( X2*X2«-Yl< Yl + Z* Z ) 41 

* ) +ATAN( X1*Y1/ABS ( Z ) /SQR T( X1*XH-Y1*Y1 +Z-*Z ) ) ) 4’ 
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APPENDIX — Continued 


DSDX< XL, X2 ,Y1 ,Y2, Z ) = Y2*AL0G( ABS ( ( X 1 + SQR T ( X1*X i+Y 2* Y2 + Z *Z ) ) / I X2 + 43 

* SQR T ( X2’i r X2+Y 2*Y2+Z*Z ) ) ) ) +Y1* ALOG ( ABS ( ( X2+SQRT ( X2*X2+Y 1*Y1+Z*Z) )/ 44 

* ( Xl+SQRT ( XI* X 1 + Y i*Y 1 i-Z’t'Z ) ) ) ) + X1 « ALOG ( ABS ( ( Y l+SQRT ( X2 *X2+Y 1 *Y 1+Z*Z 45 

* ) )/( Y2+SQPT( X2*X2+Y2*Y2+Z*Z > >*( Y2+ SORT { Xl*Xi+ Y2 *Y2 +Z *Z ) )/I Yl+SQRT 46 

* (Xl*XL+Yl*Yl +Z*Z ) ) ) l+XD*ALOG (A8S( ( Y l + SQE T ( Xl*XL +Y l *=Y 1 + Z*Z ) )/(Y2+ 47 

* SQR T ( XL *Xl +Y2*Y2+Z#Z) )*( Y2+SOP T ( X2*X2+ Y2*Y2+Z*Z ) ) / ( Y1 + SOR T ( X2*X2+ 48 

* Y1*Y1+Z*Z ) ) ») 49 

DPDY (XI , X? , Yi , Y2, Z ) = SORT (X2*X2+Y2*Y2 + Z*Z )-$*.% T (X1*X1+ Y2*Y2 + Z*Z) 50 

* -SORT (X2*X2+YI*Y1+Z*Z )+SQRT ( Xl*Xl+Y l*Y l+Z*Z ) +X1 *A LOG ( ABS ( (Xl+SQRT 51 

* (X1*XL+Y2*Y2+Z*Z ) >/( X2+SQI- T ( X2*X2+Y2* ! Y2 + Z*Z ) >*< X2+SQRT (X2*X2+Y1# 52 

* Y 1 + Z * Z ) >/(Xl+SQRT(Xl*Xi+Yi*Yl+Z*Z) ) ) )+XD*ALuG { ABS ( ( X2+SQR T ( X2* X2+ 53 

* Y2*Y2+Z*Z) 1/ (XL+SQRT (XL *XL+Y2*Y2+Z*Z ))*< XL+SQhT <XL*XL+Y1*Y1+Z*Z) ) 54 

* /{X2+S0RT ( X2*X2+Y1*Y1+Z*Z) ) ) ) 55 

D2P0XDY ( XI* X2 » Y l , Y2 * L)~ A LOG ( ABS ( (X1+SQR T (X1*X1+Y2*Y2+Z-Z) >/ (X2+ 56 

* SOR T(X2*X2+Y2*Y2+Z*Z > )*{ X2+S 0(- T (X2*X2*-Y1*Y1+Z ,? Z ) )/( Xl+SQRT (Xl^Xlt 57 

* Y 1*Y 1+Z*Z ) ) ) )+XiJ/SQR~(XL*XL+Yl*Yl + Z*Z)-XU/S,iftT(XL’i‘XL+Y2*Y2+Z*Z) 56 

OPUZt X1,X2, Yl ,Y2, Z J=Z*ALUG( ABS< ( Y2 + SORT ( X 2* X2 + Y2 *Y2 +Z* Z ) )/( Yl+SQRT 59 

* (X2*X2+Y l«Yl+Z*Z ) )*( Y1 +SQkT(X1*X1+Y1*Y1+Z*Z) ) / ( Y2 + SOR T ( XI* X 1+Y 2* 60 

* Y2 + Z*Z ) ) ) )+Xl*SGN(Z)*( ATAN< X2 + Y2/ABS (Z) /SQF T ( X2*X2+Y2*Y2 + Z*Z> )- 61 

* AT AN ( X L* Y 2/ a BS ( Z > /SORT ( X1*XL+Y2*Y2+Z*Z ) ) +A T mN ( X 1* Y 1/ ABS ( Z > /SQkT ( 62 

* X1*X1+YI*Y1 + Z*ZI )— ATANl X2* Yl / A ‘3 S ( Z ) / SQR T ( X2*X2+ Yl* Y 1 + Z* Z ) ) > + XD* 63 

* SGN { Zi*( A'AiH XL*Y2/ABS( Z) /SORT ( XL*XL+Y2*Y2+Z*Z) )-ATAN(X2*Y2/ABS(Z 64 

* ) /SQRT(X2*X2+Y2*Y2+Z*Z) ) +ATAN( X2*Y 1/ ABS l Z ) / SORT ( X2*X 2+Yl*Y 1 + Z*Z ) ) 65 

* -ATi,\j( XI_*Y1/A3S< Z J/SQRTJ XL*XL + Y L*Y l + Z*ZM) 66 

D2P0X0Z (XI, X2.Y1, Y2.Z )=SGN(Z >*< v TAN ( X 2*Y2/ ABS < Z) /SUR T ( X2*X2 +Y2*Y2+ 67 

* Z*Z> )-A T AN(Xl*Y2/A3S(Z) /SQF T ( X l *X 1+Y2*Y 2+Z* Z ) )+ AT AN( X l* Y 1 / ABS ( Z ) / 63 

* S0kT(X 1*X1+Y1*YL+Z*Z) J-ATAM X?* Y 1/AftS ( Z) /SORT ( X2#X2+Y1*Y1+Z*Z) ) ) 69 

* +XO*Z/ (XL*XL+Z*Z )*( Y2/SQRT(XL*XL+Y2*Y2+Z*Z )-Yl/SQR T ( XL*XL+Y1*Y1+ 70 

* Z*Z ) > 71 

SO ( X ) ='00 + [)R 0* ( X— X II ( 1 ) ) 72 

EL<X) = ALf.G( 1.0/S I.N(R0(X)/2.*PI ) J/4./PI 73 

DL (X ) =-OfiO*CTS(fi.v) (X)*PI/2.)/SIN(R0(X)*?I/2.)/3. 74 

PI=3. 1413925 75 

THIS R A - T OF THE PROGRAM DEFINES THE TUNNEL GEOMETRY. HERE IT IS SET UP FOR 
A SQUARE TIJNNEL OF UNIT WIDTH AND HEIGHT. EACH -WALL IS DIVIDED INTO FOUR 
STRIPS AND EACH STRIP IS DIVIDED INTO SIXTEEN ELEMENTS WHICH EXTEND FROM 
X=- -2 . 6 TO X = L . •? . 

DO l 1=1,256, 16 76 

XIHI)=-2.6 77 

XI 11 I+l )=-2.2 78 

XI l( 1+2) =-1.3 79 

XI 1( 1+3) =-1.4 80 

XI 1< I+4M-1.0 81 

X I 1 ( I +5 ) =— . 7 82 

Xlll I +6 ) =- .45 83 

XI I ( I +7 )=- .25 84 

Xll(I+3J=-.l 85 

Xll( I+9)=0.0 86 

XT 1( I +1 0 ) = • 1 87 

XI 1( 1+11 )=.25 88 

XI LC 1+12) =.45 89 

XIII I + 1 3 ) = . 7 90 

XI1U+14) = 1.0 91 

XI 1( 1 + 151=1 .4 92 

1 CON; INUE 93 
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APPENDIX - Continued 


DO 2 1=1,256 94 

XI ?C I I = X I It 1 + 1) 95 

2 CONTINUE 96 

DO 3 1=16,256,16 97 

X I 2 ( 11=1.8 98 

3 CONTINUE 99 

CO 4 1=1,64 100 

E T A1<I>=.5 101 

ETA2(I)=.5 102 

ETAl { I +64 )=-.5 103 

ETA2I 1+64)=-. 5 104 

ZETAK I +12 8 )= .5 105 

ZETA2I 1+128)=. 5 106 

ZETAK 1+192) = -. 5 107 

ZETA2 ( I +192) = -. 5 108 

4 CONTINUE 109 

DO 5 1=1, 16 110 

ZETAK I ) = .25 111 

ZETAK 1+16 1=0.0 112 

ZETAK 1 + 32 )=- .25 113 

ZETAK I +48 »=- .5 114 

5 CONTINUE 115 

CO 6 1=1,64 116 

ZETA 1 { 1+64 KZET Ai ( I ) 117 

ET All I +128 )=ZETA1 (I) 118 

ETA1( 1+192 )=ZETAl ( I ) 119 

6 COM T INUE 120' 

DO 7 1=1,128 121 

ZETA2 ( I >=Zl ; T.\l( I ) +.25 122 

FT A2 ( I + 128) = ETA1( 1+123) + . 25 123 

'< CON T INUE 124 

DO 8 1=1,256 125 

XI ( I) = ( XU ( ! ) +XI2I m /2. 12b 

ETA(I) = < ETA1 l !) +ETA2II) 1/2. 127 

ZETAII ) = ( Z t-T a 1( 1 ) +ZETA2 ( I ))/2. 126 

8 CONTINUE 129 

XTL=10.0 130 

PRINT 991 131 


THIS PART OF T“E PPUGkAM DEFINES THE WALL ChAE AC TEE i ST iCS. IN THIS CASE THE 
SIDE WALLS APE SOLID AND THE TOP AND 3 PT t nm WALLS EACH HAVE FOUR CONSTANT 
WIDTH SLOTS. THE OPEN RATIO OF T HE SLOTTED WALLS IS 6 PERCENT. 

R 00= . 06 132 

DP 0 = 0 . 0 133 

00 11 1=1, 128 134 

C 1 ( I ) =0 .0 13b 

C 2 ( I ) =0 . 0 136 

C3( I ) =1 .0 137 

C4( I ) =0 *0 136 

11 CONTINUE 139 

CO 12 1=129,256 140 

C1(!)=0.0 141 

C2( I ) = 1 .0 142 

C3( I ) =0 .0 143 

C4( T)=EL(XI(I )) 144 

12 CONTINUE 145 


17 


APPENDIX - Continued 


This PAST OF the program COMPUTES the INFLUENCE COEFFICIENTS, A(I,J) 

00 L 90 T - 1 ,266 14b 

N W I = 1 147 

1 F ( I • GT . 6 4 ) NW I = 2 148 

IF( I .G T .128 »MWI = 3 149 

I F 1 1 .GT .192)NKI =4 150 

UO 180 J = l , 25 6 151 

NWJ = 1 152 

I F ( J • G T .64 ) NW J = 2 153 

IF( J.GT . 128 )NWJ = 3 154 

IF< J.Qr .15? )NWJ=4 155 

X 1-XI ( T ) -XI 1 ( J) 156 

X2 = XI ( I )-X I2( J ) 157 

Y 1 = F T A ( I ) - t T ft 1 ( ,J ) 158 

Y2=ETA( I)-tT42( J I 159 

Z = Z£T A( I >-ZE r A( J) 160 

XU = X 1 2 ( J J — X ! I ( J ) 161 

XL=XI ( I » - X I L 162 

IFCIWJ.LT .3 >Y1=Z£TA( I )-ZE7 All J ) 163 

IF (NWJ.LT .3 )Y2=ZETAm-ZEM2 ( J ) 164 

IFC-1W J .LT .3 )Z=F TA ( I ) — E T A ( J) 165 

IFCYW I .ME .'-'WJ )90 TC 110 166 

PHi=S( X1,X2,Y1,Y2 ,1) 167 

L=US0X(X1 ,X2,Y1,Y2,Z) 168 

V=0. 0 169 

C-VCX = 0.0 170 

II = I-MUD( I, 16 ) + 1 171 

I J=J-MJD( J , lb)+l 172 

I F ( M n 0 ( I . 16) .EQ.O) IT= I-lfa + 1 173 

I F ( M 0 0 ( J, 16).FQ.0)IJ=J-16+1 174 

I F ( I .EQ.J)V=-2.U*PI*X0/2.0 175 

1F( 1 1 .EO. I J.ANO.I .G T . J ) V=-2. 0*PI*XD 176 

IF( I ,F0. J)DV0X=-2.*PI 177 

IFINW J.EQ.2)V=-V 176 

I F(NW J.EU . 3 )W=V 179 

IF (MW J.E0.4)W=-V 180 

IF (MW J.E0.2 )OVOX=-DVDX 181 

IF (NW J . EO • 3 ) OWOX=D VOX 182 

IF(NW).EQ.4)DWDX~-DVDX 183 

GO TO 130 184 

110 CONTINUE 185 

PH1=P(X1.X2, Y1,Y2,Z) 186 

U = D° 0 X ( X 1 ,X2,Y1,Y2,Z) 187 

V-CPUYl X1,X2*Y1 ,Y2,Z ) 188 

W“CP0Z(X1,X2,Y1,Y2,Z) 189 

CVDX=02PUX0Y( XI , X2, Yl , Y2 , Z ) 190 

L)WOX = 02 P OXOZ (X1,X2,Y1,Y2,Z) 191 

1 F ( M W I . G T . 2 ) GO TO 120 192 

T =V 193 

V=W 194 

W= T 195 

U T =OVOX 196 

CVUX=OaOX 197 

Dw CX = DT 198 

120 CONTINUE 199 

130 CON T IN'JE 200 

I F (NW I . Ey . 1 ) A ( I , J ) =C 1 ( I I *PHH-C2( I ) *U+C3( I )*V+C4( I) *OVDX 201 
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APPENDIX — Continued 


IFJNWI.EQ .2 ) At I , J )=C 1 ( I ) *PHI +C2( 1 ) *U-C3( I )*V-C4< I >*DVDX 202 

I F (NW I . EQ . 3 ) A ( I , J ) =C 1 ( I ) 4PHI+C2 I I >*J+-C3( 1 >*W+C4( I ) *OWDX 203 

IF(NWI .E0.4)A< I , J >=C l ( I )*PHI+C2( I ) *J-C.3( I )*W-C4( I)*OWDX 204 

180 CONTINUE 205 

190 CONTINUE 206 


the FOLLOWING STATEMENT INVERTS THE MATRIX A AND PUTS TOE INVERSE IN the 
PLACE OF THE ORIGINAL MATRIX 


CALL MATRIX<10,256,256,0, A,256,DE T > 207 


THIS PART OF THE PROGRAM COMPUTES ’Hf DISTURBANCE DUE Tu THE MODEL. HERE IT 
IS SET UP FOR A NUMBER OF VORTEX DOUBLETS LOCATED IN THE HORIZONTAL CENTER- 
PLANE OF THE TUNNEL. THE PROGRAM FIRS T READS THE NUMBER OF LIFT ELEMENTS Tu 
BE US EO AMD ThEN READS the X AND Y VALUES AND n,F wEIGnTING FACTOR FOR EACH 
ELEMENT. 


READ 94 2»L»(XA(I)»YA(I)*WT(i ),! = !, L) 208 

DO 200 1=1,256 209 

Bill =0.0 210 

200 CONTINUE 211 

DO 299 K= 1 » l 212 

XPP=XA(K) 213 

YPP = Y A ( K ) 214 

CO 210 1=1,64 215 

X=XI(I)-XPP 216 

Y = E t A( I ) - YP P 217 

Z=ZETA(I) 218 

PHI=Z/(Y*Y+Z^Z)’M 1.0+X/SJRT< X*X + Y*Y+Z*Z ) ) /4.0/PI 219 

L = Z/4. / P I / ( X*X+Y*Y+Z*Z)**1.5 220 

V=-( 2.* Y*Z + ( 2 . *X* >3^Y*Z *3 ,»x*Y* :• 3* Z «- 3 .*X~ Y*Z**3) / ( X*X+Y*Y+Z *Z ) <'*1 . 221 

* 5J/4./PI / (Y*Y+Z*Z)* : *2 222 

DVDX=-. 75/P !*Y*Z/ ( X*X + Y*Y+Z*Z ><'*2.5 223 

e ( I ) = B( I ) +( Cl ( I )*PHH-C2 m *LM-C3( I )* V«-C4C)*0VDX)*WT (K ) 2 24 

210 CONTINUE 225 

DO 220 1=65,1.20 226 

X = X I ( I )-XPP 227 

Y=E t a( I ) — YPP 228 

Z = ZET A ( I ) 229 

PHI*Z/CY*Y + Z*Z)*(1.0-HX/SQRT(X*X + Y*Y+Z*Z))/4.0/PI 230 

U-Z/4./PI/( X*X+Y*Y+Z*Z)**i.5 231 

V=-(2.*Y«Z+(2. *X* *3*Y*Z+3 . *X*Y** 3*Z*2.4X4Y*Z443)/{X4X+Y*Y+Z4Z)**l. 2 32 

* 5)/4./PI/CY*Y«-Z*Z)**2 233 

D VDX=— . 75 / P 1 ■•-Y*Z/ ( X*X+Y*Y + Z* Z )**2. 5 234 

E(I)=B(I) + (C1CI ) * P H i +C 2 ( I >*U-C3 ( I )*V-C4( I )*DVDX)*WT (K ) 235 

220 CONTINUE 236 

CO 230 1=129, 192 237 

X=XI(I)— XP» 238 

Y = ET A ( I ) -YPP 239 

Z= ZET A ( I ) 240 

PH I = Z / ( Y*Y+Z*Z)*( 1 .0+X/SORTI X*X+Y*Y+Z*Z ) J/4.0/PI 241 

U=Z/4./P I/( X*X+Y*Y+Z*Z) **1 .5 242 

W=(Y*Y— Z*Z+(X** 3* Y4Y+X-Y* •=4-X*43*Z*Z-X*Y*Y*Z*Z-2 .*X*Z**4 ) /< X*X+Y*Y 243 

* +Z*Z)**1.5)/4./PI/(Y*Y+Z*Z)**2 244 

DWDX = ( X* X *- Y * Y— 2 . * Z *Z ) / { X*X+Y*Y*Z *Z ) **2 . 5/4 . /P I 245 
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APPENDIX — Continued 


e(I)=6(I)+(ClU )*PHI+C2 t I )*U+C3( I ) *W+C4 ( I ) *DWDX) *WT(K) 246 

230 CONTINUE 247 

DO 240 1=193,256 248 

X=XI(I)-XPP 249 

Y = £ T A ( I ) — Y PP 250 

Z=ZETA( I ) . 251 

PH I =Z/ ( Y*Y+Z*Z ) # ( l .O+X/SQRT ( X*X*Y*Y«-Z*Z ) ) /4.0/PI 252 

U=Z/4./PI/( X*X+Y*Y *Z'*Z )##1 .5 253 

W= ( Y *Y— Z*Z+ ( x**3* Y*Y-*-X*Y **4— X#< ! 3*Z*Z-X !!s Y*Y J i £ Z-'i £ Z-2 ,*X*Z** 4) /( X*X<-Y*Y 254 
* +Z *Z)**1.5)/4./PI/( Y*Y+Z*Z)**2 255 

DWDX = ( X*X+Y*Y-?.*Z*Z )/( X*X+Y*Y+ Z*Z >** 2. 5/4. /P I 256 

6(1 ) = B( I ) + ( C 1 ( I )*PHI+C2( I )*U-C3( I ) * W-C4 ( I )*L)NOX)*WT ( K ) 257 

240 CONTINUE 258 

259 CONTINUE 2 59 


THIS PART HE T HE PH OOP AM COMPUTES THE SOURCE STRENGTH SLOPES WHICH SATISFY 
THE BOUNDARY CONDITIONS. 


DO 300 1=1,256 260 

SIGMA! I ) = 0 . 0 261 

300 CONTINUE 262 

DO 302 1=1,256 263 

DO 301 J = 1 , 2 5 6 264 

SI GMA( I)=S IGMA( T )-A(I , J )*B( J > 265 

301 CONTINUE 266 

302 CONTINUE 267 


THIS P A 0 T UP THE PROGRAM COMPUTES THE LIFT INTEkFERENCE FACTOR. IT READS 
THE X, Y, AND Z VALUES AT WHICH THE INTERFERENCE IS TO 8E COMPUTED. 


CONTINUE 

268 

READ 993 , XC , YC, ZC 

269. 

IFIEOF, 5)990, 401 

270 

CONTINUE 

271 

CELTA1=0.0 

272 

CELT A 2 = 0 . 0 

273 

CELTA3=0.0 

274 

DELTA4=0.C 

275 

DO 410 J = 1 , 64 

276' 

X1 = XC-XI 1 ( J ) 

277 

X2 = XC--XI2( J) 

278j 

Y=YC-ETA( J ) 

279 

Z1 = ZC-ZET41 { J ) 

280 

Z2 = ZC-ZETA2< J ) 

281 

X D = X 1 2 ( J) — XII (JJ 

282 

XL =XC-X I L 

283 

DP D Y (XI ,X2,Z1,Z2,Y) 

284 

DELT Al = CEL T Al +W* SIGMA ( J)/2. 

285 

CL NT I NUE 

286 

DO 420 J= 65 , 1 2 8 

287 

X1 = XC-XI 1( J ) 

288 

X2 =XC— X I 2 ( J ) 

289 

Y = YC— ET A ( J ) 

290 

Z 1 =ZC— Z ET A 1 ( J I 

291 

Z2 = ZC-ZETA2< J J 

292 
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APPENDIX — Concluded 


Xl)=X I 2( J >— X T 1 ( J) 293 

XL =XC-X I L 294 

V. = bP0Y(Xl,X2,Zl ,Z2,Y) 293 

CELT42=D£LTA2*W#SIGMA( J)/2. 296 

A 20 CONTINUE 297 

DO 430 J - 1 2 9 1 192 298 

X1=XC-XI 1 ( J ) 299 

X2 = XC-X 1 2 ( J ) 300 

Y1=YC— ETA1( J) 301 

Y2 = YC-ETA2U> 302 

Z=ZC-ZETA( J ) 303 

XD = X I 2 ( J ) -X I 1 ( J ) 304 

XL =XC-X I L 303 

W=DPDZ(X1,X2,Y1,Y2,Z) 306 

DELTA3 = DbLTA3+W'!'SIG.MA(J)/2. 307 

430 CONTINUE 308 

DO 440 J •= 1 9 3 * 256 309 

Xl=XC-XI 1( J) 310 

X2 = XC— X I 2 ( J ) 311 

Y 1 = YC- ET A 1 ( J ) 312 

Y2 = YC — ETA2I J) 313 

Z=ZC-ZETA ( J) 314 

XD=X I 2 { J )-XI 1 ( J ) 313 

XL=XC-XIL 316 

Vv=CPDZ< XI, X2» Y1 * Y2 * Z ) 317 

CELT A4=DELT A4+W*S IGMA( J ) / 2 . 318 

440 CONTINUE 319 

CELTA = OELTA1+OELTA2 + L)ELT A3+DELT a 4 320 

PRINT 993, XL, YL,ZC., DELTA 321 

GO T 0 400 322 

990 STOP 323 

991 FORMAT ( 39H1 X Y Z DELTA) 324 

992 FORMAT! I2,/3E10.6) 325 

993 FCPMAT ( 4,F1 0 .6 ) 326 

END 327 
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